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Abstract 

Given an RNA sequence a, consider the network G = {V,E), where the set V of nodes consists of all 
secondary structures of a, and whose edge set E consists of all edges connecting two secondary structures 
whose base pair distance is 1. Define the network connectivity, or expected network degree, as the average 
number of edges incident to vertices of G. Using algebraic combinatorial methods, we prove that the 
asymptotic connectivity of length n homopolymer sequences is 0.473418 • n. This raises the question of 
what other network properties are characteristic of the network of RNA secondary structures. Programs 
in Python, C and Mathematica are available at the web site http: //bioinformatics .be.edu/clotelab/ 
RNAexpNumNbors 


1 Introduction 

In [24], Stein and Waterman used generating function theory to determine the asymptotic number of RNA 
secondary structures. Since that pioneering paper, a number of results concerning RNA secondary structure 
asymptotics have appeared, including the incomplete list [ISl [HI UHl 131 HZl 13 Ell EH 1311313113 E3 ESI IH] ■ 

In contrast to the previous papers, here we consider network properties of the ensemble of RNA secondary 
structures, following the seminal work of Wuchty ESI, who showed by exhaustive enumeration of the low 
energy secondary structures of E. coli phe-tRNA, that the corresponding network architecture had small- 
world properties |5S|. Small-world networks appear to abound in biology, providing a kind of robustness 
necessary for the molecular processes of life, as seen in networks of neural connections of C. elegans [29] . 
gene co-expression in S. cerevisiae [23, metabolic pathways E3I2I], intermediate conformations in tertiary 
folding kinetics for the protein villin [T], etc. 

In this paper, we use algebraic combinatorial methods, and in particular the Flajolet-Odlyzko the¬ 
orem |7] to prove that the asymptotic expected network connectivity of RNA secondary structures is 
0.4734176431521986 • n. Following [351 E3], stickiness is defined to be the probability p that any two po¬ 
sitions can pair. For the simplicity of argument, in the homopolymer model, we take stickiness p to be 1; 
however, minor changes in our C dynamic programming algorithm and in our Mathematica code permit the 
computation of asymptotic expected connectivity for arbitrary stickiness p. 


2 Preliminaries 

A secondary structure for a given RNA nucleotide sequence a = ai,..., a„ is a set s of base pairs {i,j), 
such that (i) if {i,j) G s then a, a^ form either a Watson-Crick (AU,UA,CG,GC) or wobble (GU,UG) base 
pair, (ii) if (j, j) S s then j — i>9 = 3{a steric constraint requiring that there be at least 0 = 3 unpaired 
bases between any two paired bases), (Hi) if {i,j) € s then for all / ^ j and i' ^ i, {i',j) ^ s and (i,j') ^ s 
(nonexistence of base triples), (iv) if {i,j) € s and {k,£) G s, then it is not the case that i < k < j < i 
(nonexistence of pseudoknots). For the purposes of this paper, following Stein and Waterman [33], we 
consider the homopolymer model of RNA, in which condition (i) is dropped, so that any base can pair with 
any other base. 

Suppose that a = ai,..., a„ is an RNA sequence. If s is a secondary structure of a, then let denote 
the number of secondary structures of a that can be obtained from s by the removal or addition of a single 
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base pair; i.e. those structures having base pair distance from s of 1. Define the expected number of neighbors 
(Ns) by 

m = I (1) 

where Q = total number of neighbors of all secondary structures s of a, and Z de¬ 

notes the total number of secondary structures of a. Note that Z corresponds to the partition funetion 
Yfs exp{—E{s)/RT) if the energy of every structure is 0. 

In [1] we described three algorithms to compute the expected number of neighbors, or network connec¬ 
tivity, (Ns) = exp{-E{s)/RT) ^ where Z is the partition function exp(—ill(s)/i?T) with respect to 

energy model A (each structure s has energy E{s) = 0), model B (each structure s has Nussinov energy 
E{s) equal to —1 times the number |s| of base pairs in s), and model C (each structure s has energy E{s) 
given by the Turner energy model naisi]). 

Below, we follow reference [4] in deriving the recurrence relations for Q and Z for model A, corresponding 
to equation equation Q. For 1 < i < j < n, define the subsequence a[f, j] = a^,..., a^, and define S'5'(a[f, j]) 
to be the collection of secondary structures of sd(i,j]. Define 

Q^,3 = E (2) 

sGSS(a[i j]) 

Similarly, let Zi^j = SsGSS(a[i j]) i-®- denotes the number of secondary structures of a[z, j]. 

Base Case: For j - i e {0,1, 2, 3}, = 0 and Z^j = 1. 

Inductive Case: Let BP{i,j, a) be a boolean function, taking the value 1 if positions i,j can form a base 
pair for sequence a, and otherwise taking the value 0. Assume that j — i > S. 

Subcase A: Consider all secondary structures s G a[f, j], for which j is unpaired. For each structure s in 
this subcase, the number Ng of neighbors of s is constituted from the number of structures obtained from 
s by removal of a single base pair, together with the number of structures obtained from s by addition 
of a single base pair. If the base pair added does not involve terminal position j, then total contribution 
to J2s£SS(sL[i j]) It remains to count the contribution due to neighbors t of s, obtained from 

s G S'S'(a[z, j]) by adding the base pair {k,j). This contribution is given by X]i=i BP{k, j, a.)-Zi^k-i'Zk+i,j-i, 
where Zi^i-i is defined to be 1. Thus the total contribution to Qij from this subcase is 

t-4 

Qi,j—i T ^ ^ BP(^k, j, a) • Zi j„—i ' Z}^^i j—i. 

k—i 

Subcase B: Consider all secondary structures s G a[z,j] that contain the base pair (fc, j) for some k G 
{i,... ,j — 4}. For secondary structure s in this subcase, the number Ng of neighbors of s is constituted 
from the number of structures obtained by removing base pair {k,j) together with a contribution obtained 
by adding/removing a single base pair either to the region [i, fc — 1] or to the region [k -|- 1, j — 1]. Setting 
Qi^i-i to be 0, these contributions are given by 

t-4 

k—i 

In the current subcase, the contribution to Zij is BP{k,j,a) ■ Zi^k-i ■ Zk+ij-i- 

Finally, taking the contributions from both subcases together, it follows that 

I-4 

Qi,j — Qi,j—i F ^ ( BPi^k, j, a) • [2 • -t- Qi^k—i * i F Z-ij^—i • (3) 

k—i 

i-4 

Zi,j = Zjj-i + BP{k,j, a) ■ Zj^k-i • (4) 

k—i 
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Figure 1: (Left) All possible secondary structures of the 7-mer homopolymer, where position i can pair with 
position j provided only that 1 < i < j < 7 and j i > 4. (Right) Graph representation of neighborhood 
network, where nodes a,b,c,d,e,f,g,h respectively represent the 8 secondary structures in the list. The number 
of neighbors of each secondary structure is indicated to its right. The expected number of neighbors for the 
7-mer is thus 2)/8 = 16/8 = 2. 


It follows that the expected number ( Ns ) of neighbors Ns of structures s of a is Note that the recursion 

for Zij is well-known and due originally to Waterman m- 

To provide concrete intuition for the problem we consider, in Figure 1, we present the list of secondary 
structures for a homopolymer of length 7, depicted as a network having expected connectivity of 2. In the 
left panel of Figure 2, we present a histogram for the network connectivity (graph degree or number of 
neighbors), by analyzing an exhaustively produced list of all 106,633 structures of the 20-mer homopolymer. 
In the right panel of Figure 2, we present a plot of the normalized expected number of neighbors of for 
homopolymers of length 1 to 1000 nt, obtained by dividing the expected number of neighbors by sequence 
length. Clearly there appears to be an asymptotic value for the length-normalized expected connectivity, 
suggesting that it may be possible to formally prove the existence of this asymptotic value, a task to which 
the remainder of the paper is dedicated. 

3 Methods 

3.1 Recurrence relation for partition fnnction 

In order to determine asymptotic network connectivity, we now provide similar recursions to those of (§ 
and Q for the homopolymer model of RNA, where any base (position) i can pair with any other position 
(base) i, provide only that l<i-|-0-|-l<j<n. Following common convention due to steric constraints, 
we take 6 to be 3. These recursions are the basis of the dynamic programming code we implemented, used 
to produce Figure 

We define Zq = 1, in order to simplify the recurrence relation for Z„, defined to be the number of 
secondary structures for a homopolymer of length n, or equivalently, the partition funetion for the energy 
model that assigns an energy of 0 to every structure. Moreover, since the empty structure is the only 
structure for sequences of length 1, 2, 3,4 = 0 -|- 1, we define = 1 for 0 < n < 4. Secondary structures for 
a homopolymer of length n > 4 can be partitioned into two classes: (1) n is unpaired, (2) there is a base 
pair (fc, n) for some l<fc<n — 4 = n — 0 — 1. Thus we have 

_ri if0<n<4 = 6»-hl 

\ Z„_i -f Zk ■ ^n-2-fc if n > 5 = 0 + 2 

which is the homopolymer analogue of equation Q. To employ generating function theory, we require a 
single formula for rather than a definition by cases - see p. 66 of [9]. This is easily achieved by (i) 
adding the indicator function I[n = 0], defined to equal 1 if n = 0, and otherwise 0, and (ii) adding and 
subtracting the same terms to ensure that k ranges from 0 to n — 2, rather than n — 9 — 2 = n — 5. Thus 
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Plot of normalized expNumNbors as function of sequence length 



Figure 2: (Left) Histogram of the number of neighbors for all 106,633 secondary structures of the 20-mer 
homopolymer. The mean is 8.336, the standard deviation is 4.769, the maximum is 136, and minimum is 
3. (Right) Plot of the normalized expected number of neighbors of for homopolymers of length 1 to 1000 
nt, obtained by dividing the expected number of neighbors by sequence length. Apparent asymptotic value 
seems to be « 0.4724. The main result of this paper is the proof that the asymptotic value is in fact 
0.4734176431521986. 

equation ([^ is equivalent to ^ defined by: 


Zn — Z, 


n — 1 


Zk ■ Zn-2-k + — 0] — ^n-1 ■ — Zn_^ • Zi — ■ Z 2 

n —2 

— Zji—i + ^ ^ Zk ■ Zfi—2—k Zji—\ Zf^—^ Zfi—4^. 


fc =0 


Let z = multiplying equation (|^ by a;”, summing from n = 0 to 00 , we obtain 

00 00 00 n —2 00 

Y,Zn-X^ = Y, Zn -1 • X" + ^ ^ Zfc • Z„_ 2 -fe • x" - ^ (Z „_2 + Z„_3 + Z„_4 - I[n = 0]) • x" 

n—Q n—0 n—0 k—0 n—0 

which yields 

I 2 2 2 3 4,1 

z = xz + xz— zx— zx — zx +1 
Solving the quadratic equation Q for z, we determine that 


z = 

Only the first solution 


1 — X + x^ + x^ + x"* ± \/l — 2x — x^ + x^ + 3x® + 2x'^ + : 


2x2 


zl = 


1 — X + x^ + x^ + x^ — Vl — 2x — x2 + x^ + 3x® + 2x^ + : 

2^2 


( 6 ) 


(7) 


( 8 ) 


has the property that the coefficients of its Taylor expansion correspond to the values of as determined 
by dynamic programming (see program at web site). In particular, using Mathematica, we obtain 

zl = 1 + X + x^ + x^ + x^ + 2x^ + 4x® + 8x’' + 16x® + 32x® + fibx^^ + 133x^1 + 274x^2 + 568x^3 + 

1184x^'‘ + 2481x^'^ + 5223x^® + 11042x^^ + 23434x^® + 49908x^® + 106633x^° + 0{xf^ 
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which can be compared with the values determined by our dynamic programming implementation: 


n 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

Zn 

1 

1 

1 

1 

1 

2 

4 

8 

16 

32 

65 

133 

274 

568 

1184 

2481 

5223 

11042 


3.2 Recurrence relation for the number of neighbors Qn 

Let Qn denote the total number Ng of nearest neighbors of all secondary structures s of a homopolymer of 
length n. The homopolymer analogue of equation (|^ is as follows: 


Qn 


0 


if n = 0,1, 2, 3,4 


Qn-l+^k=0 ‘^^kZn-2-k + QkZn-2-k + ZkQn-2-k else. 


(9) 


If we assume that Zn = 0 = Qn for n < 0, then it follows from equations and ([^, that we can express 
Qn by the formula 


n—5 


Qn — Qn—1 ^ ^ ‘^^k^n—2—k Qk^n—2—k ^kQn—2—k- 


( 10 ) 




Next, we add and subtract the same terms in order to ensure that the upper bound in the previous summation 
is n — 2. This yields 


n-2 


Qn — Qn—1 T E ‘^ZkZn-2-k + QkZn-2-k + ZkQn-2-k 


k—Q 


~2 {Zn-2 • Zq + Z „_3 • Z\ + Z „_4 • Zq) — {Qn-2 ' Zq + Qn-3 ' Z\ + Qn-i ’ Z 2 ) 
— {Zn-2 • Qo + Zn-3 ' Ql + -^71-4 ’ Q 2 ) 


which simplifies to 


n-2 


Qn — Qn-1 + ‘^ZkZn-2-k + QkZn-2-k + ZkQn-2-f. 


k=0 


~2 {Zn-2 + Zn-3 + Zn-i) — {Qn-2 + Qn-3 + Qn-i) ■ 


Multiply each term of equation (12) by and summing from n = 0 to cx), to obtain 

oo oo oo / n—2 ^ 


oo / n—2 


E! Qn • X — Qn-1 • X + 2ZkZn-2-k | + E] I E] Qk^n — 2—k j T 

n—0 n—0 n—0 \k—0 / n—0 \k—0 / 

oo /n—2 \ oo oo 

EE ZkQn—2 — k -E 2 {Zn-2 + Zn-3 + Zn-i) — {Qn-2 + Qn-3 + Qn- 4 ) • 


n—0 \k—0 


n—0 


n—0 


Let q = Qn^;" and z = ZnX'^- Then from (13) we have 

q = xq + 2x^z^ + x'^qz + x^zq — 2x^z — 2x^z — 2x‘^z — qx^ — qx^ — qx'^. 

From equations ([^ and (14), we have 

• z'^x'^ = z — zx + zx^ + zx^ + zx^ — 1 

• q = xq + 2a;^z^ + 2x'^qz — 2x^z — 2x^z — 2x'^z — qx'^ — qx^ — qx^ 

from which we eliminate the variable z to obtain the following quadratic equation in variable g, 

4a:® = q^x^ {l — 2x — x^ + x"^ -\- 3a;® + 2a;^ + a;®) + 

g (2 - 6a: + 2a;^ + 2x® + 2a;^ - 2a;® + 6a:® - 2a:^ - 2a;® - 2a:®) . 


( 11 ) 


( 12 ) 


(13) 


(14) 


(15) 
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Solving for g, we determine that only the solution 


q2 = (^—l + ^x — x^ — x^ — x‘^ + x^ — ix^+x^ + x^ + x^+ (16) 

y/ (l — 6x + lla;^ — Ax^ — Zx'^ — 6a;® + 15a;® — 16a;^ + a;® + 2a;® + 9a:^®— 

8x“ + Ix^^ + 6x^3 + bx^'^ + 3x^® + 2x^^ + x^®)) / (x^ - 2x® - x^ + x® + 3x® + 2x® + x^®) 

is possible, since its Taylor expansion is 

q2 = 2x® + 6x® + 16x^ + 40x® + 96x® + 228xi° + 532x“ + 1230xi® + 2826x^® + 6464xi'‘ + 14742x^® + 
33546x1® + 76216x11' + 172968xi® + 392228xi® + 888932x®° + Oix)^\ 


These values agree with those from the following table that were obtained by our dynamic programming 
implementation. 


n 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

Qn 

0 

0 

0 

0 

0 

2 

6 

16 

40 

96 

228 

532 

1230 

2826 

6464 

14742 

33546 

76216 


3.3 Flajolet-Odlyzko theorem 


we use the following 
resp. Qn] in the Taylor 


Having determined the formulas for zl [resp. < 72 ] in equation ([^ [resp. equation (16 
theorem of Flajolet and Odlyzko [7] to determine the asymptotic coefficients of Z, 
expansion of the generating function zl [resp. q2]. 

Following standard convention, let [x"']/(x) denote the nth coefficient in the Taylor expansion of /(x 
The following theorem is stated as Corollary 2, part (i) of [7] on page 224. 


Theorem 1 (Flajolet and Odlyzko) Assume that /(x) has a singularity at x = p > 0, is analytic in the 
rest of the region A\l, depicted in and that as x ^ p in A, 


f{x)^K{i-x/pr. 


(17) 


Then, as n ^ oo, if a ^ 0, 1 , 2,..., 

/„ = [x"]/(x) ^ • n-“-V- 

r(-a) 

where F denotes the Gamma function. 


we determine the asymptotic value of [x^Jzl = Z^- The ratio of these values then yields the asymptotic 
expected number of neighbors, or network connectivity for the homopolymer model of RNA. 


In Section 3.4 we determine the asymptotic value of [x'^]q2 = and in the following Section 3.5 


3.4 Asymptotic number of neighbors 


Let P denote the polynomial under the radical of equation (16), i.e. 


P = 1 — 6 x + llx® — 4x® — 3x‘^ — 6 x® + 15x® — 16x^ + x® + 2x® + 

9 x 1 ° _ g^ll 7^12 g^l 3 53,14 33 .I 6 2 x 11 ’ + xi®. 


(18) 

(19) 


There are 4 real roots and 14 imaginary roots of P; however, the root of the smallest modulus (absolute 
value) is the real root 


p = 0.436911127214519 « 0.436911. (20) 

Now q2 can be expressed in the form q2 = G + H, where G = and H = , and where 

Gnum = (—1 + 3x — X® — X® — x"i + X® — 3x® + xi” + X® + X®) (21) 

Gdenom = x® (l — 2x — x® + x^ + 3x® + 2x1' + x®) . 
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Figure 3: The shaded region A where, except at a: = p, the generating function f{x) must be analytic. Here 
p = 1. Figure taken from Lorenz et al. HZl. 


One notes that p is a root of both Gnum and Gdenom, so their ratio is well-defined; as well, clearly 0 is 
a root of Gdenom. It follows that both 0 and p are singularities of the function q2 (recall that in complex 
analysis, the square root of 0 is a singularity). 

For asymptotics, the term can be neglected, since lima;_>p cde^m ~ —2.963617606602476; thus 

it follows that [x'^\q2 = [x^] ■ However, the Flajolet-Odlyzko theorem cannot be applied to the 

function / = ; since p is not the dominant singularity, as 0 < p. To address this issue, we define 

Gnuml = Gnum and Gdenoml = Gdenom/x'^, hence 

Gnuml = (^—l + ix — x^ — x^ — x'^ + x^ — Zx^+x'^ + x^ + x'^) (22) 

Gdenoml = {l — 2x — x^ + x*" + 3x® -I- 2x'^ + x^) . 

Now we can apply Theoremto the function fl = x^-f= , for which p is the dominant singularity. 

First, we factor (1 — x/p) out from P. 


- -- = 1 - 3.71121a; -h 2.50581a;2 -h 1.73529a;3 0.971726a;^ - 3.77592a;® -h 6.3577a;® - 

1 - x/p 

1.44854a;^ - 2.3154a;® - 3.29948a;® -f 1.44816x^° - 4.68545a;^^ - 3.72404a;^® - 
2.52356a;^® - 0.775919a;®'‘ - 1.77592a;^® - 1.06471x^® - 0.436911x®^. 


Second, we factor (1 


Gdenoml 
1 - x/p 


— x/p) out from Gdenoml. 

= 1 + 0.288795X - 0.339007x® - 0.775919x® - 0.775919x‘‘ - 1.77592x® - 
1.06471X® - 0.436911x^. 


(23) 

(24) 

(25) 


(26) 


It follows from (23) and (26) that 


Vp 

Gdenoml 
^P/il-x/p) 
Gdenoml / {1 — x/p) ^ 


(1 - x/p) ■ yJP/{l - x/p) 
Gdenoml/[1 — x/p) 

0.11422693623949792. 


(27) 

(28) 


Gd^Zmi = 0.11422693623949792 • (1 - x/p)-^/^. If we define K = 0.11422693623949792, then it 
follows that fl = ~ by applying Theorem for a = —1/2, we have the 

following asymptotic result. 
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Lemma 1 The asymptotic value of the nth coefficient [x^] (^x^q2) in the Taylor expansion 

0.06444564758689844-2.2887949921884863^ 

\/n 


IS 


Proof. We have 


[x-](x^q2) = [x"](/l) 


K 


.p-" 


r(-a) 

0.11422693623949792 ^ r, 

■---• 0.436911127214519""' 

r(i/2) 

0.06444564758689844 • 2.2887949921884863"^ 


( 29 ) 


3.5 Asymptotic number of structures 

We now proceed similarly to determine the asymptotic value of the Taylor coefficients of x^ ■ 

Now let P denote the polynomial under the radical of equation Q, i.e. 

P = 1 — 2x — + 3a:® + 2x^ + a;®. (30) 

There are 2 real roots and 6 imaginary roots of P; however, the root of the smallest modulus (absolute value) 
is the real root 


p = 0.436911127214519 « 0.436911, 


(31) 


identical the the value from equation (|20). 

Vp 


G = 


Gnum 

Gdenom 


and H = 


Gdenom ’ 


and where 


Then zl can be expressed in the form zl = G + H, where 


Gnum = (l — X + + X® + x^) (32) 

Gdenom = 2x^. 


Note that p is not a root of either Gnum or Gdenom, so their ratio is well-defined; as well, clearly 0 is a root 
of Gdenom. It follows that both 0 and p are singularities of the function zl (recall that in complex analysis, 
the square root of 0 is a singularity). 

For asymptotics, the term can be neglected, since lima;_,.p = 2.2887949921884205; thus 

it follows that [x"]zl = [x"] ■ However, the Flajolet-Odlyzko theorem cannot be applied to the 

function / = since 0, rather than p, is the dominant singularity. To address this issue, we define 

Gnuml = Gnum and Gdenoml = Gdenom/x"^, hence 

Gnuml = (l — X + x^ + X® + x^) (33) 

Gdenoml = 2. 

Now we can apply Theoremto the function fl = x^-f= Q,Tilnni ’ which p is the dominant singularity. 
First, we factor (1 — x/p) out from P. 


P 


1 - x/p 


= 1 -f 0.288795X - 0.339007x2 _ o.7759i9a;3 _ 0.775919x‘‘ - 1.77592x® - 1.06471x® - 0.43691|fM) 


It follows from (34) that 


Vp 


(1 - x/p)l/2 . ^P/(l-x/p) 


Gdenoml 

VP! a -x/p) 


Gdenoml 


ip) = 


VP HI-x/p) 


(p) = 0.4825630725501931 


(35) 

(36) 


Gdenoml '' ' 2 

and soVf = 0.4825630725501931 • (1 - x/p)i/2. if we define K = 0.4825630725501931, then it follows that 
fl = K{1 — x/p)^/2^ and by applying Theorem for a = +1/2, we have the following asymptotic 


result. 























Lemma 2 The asymptotic value of the nth coefficient [x'^] (x^zl) in the Taylor expansion of x'^ 


IS 


0.13612852946880957-2.2887949921884863" 


ri3/2 


Proof. We have 


(x^z 


1 ) = 


(/I) 


K 


n-a) 

0.13612852946880957 




• 2.2887949921884863” 


r(-i/2) 

0.13612852946880957 • 2.2887949921884863” 

j-i3/2 


(37) 


□ We have now established the asymptotic expected network connectivity. 


Theorem 2 The asymptotic value of the expected number of neighbors is 0.4734176431521986 • n; i.e. the 
asymptotic value ofJ2^=o ** 0.4734176431521986. 

Proof. By Lemmas [T] and we have 

[ai”]g2 [a;”](a;2 • q2) 

[a;"]zl [x'^]{x‘^ ■ z\) 

(0.06444564758689844 • 2.2887949921884863” • r-^/^ 

0.13612852946880957 • 2.2887949921884863” • n-3/2 
= 0.4734176431521986 • n. 
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